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Abstract In this article, a one-dimensional and a two- 
dimensional approach to the evaluation of local diffusion 
coefficients for Norway spruce sapwood from measured 
moisture content (MC) values are presented. A studied 
wood sample was dried from the initial green condition to 
about 15% mean MC, but here only the diffusive part of the 
drying process between approximately 25% and 15% mean 
MC was treated. Measured local MC values were based on 
nondestructive X-ray computed tomography data. Finite 
element calculations were performed with two alternative 
diffusion coefficients to test the appropriateness of the dif¬ 
fusion coefficients that were evaluated from the measured 
MC values. The evaluated diffusion coefficients show inter¬ 
esting dependence on MC and distance from the evapora¬ 
tion surface. The advantage of using the methods presented 
is that the diffusion coefficient is calculated on a local level 
without having to define a function for the diffusion 
coefficient’s dependency on other parameters. 

Key words Wood • FEM • Drying • Diffusion coefficient • 
Computed tomography 


Introduction 

When physically describing the drying behavior of wood, 
the drying process can be divided into the capillary part and 
the diffusion part. In the capillary part, the moisture content 
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(mass of water/mass of dry wood; MC) is high and there is 
free water present in the voids of the wood fibers. When the 
MC is lower and there is moisture only bound in the cell 
walls of the fibers, then the moisture flux is driven by diffu¬ 
sion. Internal stresses are induced in the diffusion part of 
the drying process due to anisotropic shrinkage, which may 
cause checking and drying distortions that reduce the qual¬ 
ity of the timber. Large MC gradients during the diffusive 
process give a fast drying process but cause large stresses. It 
is of importance to understand the drying behavior of wood 
in order to avoid quality degradation due to drying. One 
way to do this is through simulation. 

Simulations of the wood drying process using three- 
dimensional (3D) finite element method (FEM) can pro¬ 
vide detailed and realistic information about the local MC, 
local stress, and global deformation history. 3D FEM calcu¬ 
lations require, among other material data, diffusion coeffi¬ 
cients that are valid locally throughout the material. The 
objective of this work is to determine local diffusion coeffi¬ 
cients using nondestructive measurements and numerical 
methods. 

Alternative approaches to the evaluation of local diffu¬ 
sion coefficients for Norway spruce sapwood are presented 
based on experiment. A clearwood sample (Fig. 1) was 
dried from the initial green condition to about 15% mean 
MC (mean value for the wood sample in question), but here 
only the diffusive part of the drying process between ap¬ 
proximately 25% and 15% mean MC was treated. The mea¬ 
sured local MC values were based on nondestructive X-ray 
computed tomography (CT) data. 

Several authors (Hukka, 1 Hukka and Oksanen, 2 Liu 
et al., 3 Rosenkilde and Arfvidsson 4 ) have found that the 
diffusion coefficient for a certain wood sample is not con¬ 
stant, but is dependent on MC in addition to the depen¬ 
dence on temperature. The diffusion coefficients obtained 
are global mean values for a wood sample of a certain size. 

The CT method used here to measure local interior two- 
dimensional (2D) densities and MCs of the wood sample is 
described by Lindgren, 5 Danvind and Moren, 6 and Wiberg. 7 
The advantage of using the methods presented in this article 
is that the diffusion coefficient is calculated on a local level 
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are attracted to wood molecules and not to a specific vol¬ 
ume. Mass conservation and Eq. 1 give 
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Fig. 1. Wood sample. Asterisks and circles denote center points of cells 


without having to define a function for the diffusion 
coefficient’s dependencies on other parameters such as MC. 
The possibility of using CT data to determine local proper¬ 
ties in wood other than moisture properties, such as spiral 
grain angles, is demonstrated by Ekevad 8 and Sepulveda 
et al. 9 
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which is a modified Fick’s second law. The general bound¬ 
ary conditions to be used with Eqs. 1 and 2 are specified 
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(convective conditions or mixed conditions) on all or parts 
of the boundaries. Here the convective condition is used for 
the evaporation surface and the natural condition g = 0 for 
all the other surfaces. The initial condition is u- u 0 (x , t 0 ) at 
the starting time t Q . ft is the mass transfer coefficient for the 
moisture vaporization into the surrounding air on the 
boundary surface, and u M is the equilibrium MC for wood 
under ambient air conditions. 

The equations for 2D mass flux in the x and y directions 
are a modified Fick’s first law for an orthotropic material, 
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where 


Theory 


Isothermal conditions are considered, and a Cartesian co¬ 
ordinate system and a referential (Lagrangian) viewpoint 
are adopted in order to deal with the shrinkage of the 
material. Thus, the coordinates x, y, and z denote the coor¬ 
dinates of a material point in the green condition, and all 
lengths, areas, and volumes are values in this green state 
and remain constant during the drying process. In the first 
approach to determining the diffusion coefficient D, the 
moisture flux is assumed to be one dimensional (ID) in the 
radial direction (x) and a modified Fick’s first law is 
stated as 
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is the symmetric diffusion coefficient matrix. Here it is 
assumed that radial and tangential directions coincide 
with the x and y directions, respectively (Fig. 1). Mass con¬ 
servation and Eq. 4 give 
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where g = g(x, t ) is the mass flux in the positive x direction 
per unit area at position x at time t, u = u(x , t) is the MC at 
x and t. u is used as the driving potential in the modified 
Fick’s first law (Eq. 1) instead of using the moisture concen¬ 
tration w-upQ, where p 0 = p 0 (x) is the basic density of wood 
(dry mass of wood per green volume). It was found experi¬ 
mentally that u is a better way to express the amount of 
moisture in our case where p 0 varied in space (Fig. 2). This 
was confirmed in our experiments, because it was found 
that a gradient dw/dx / 0 could exist without creating a mass 
flux, due to local variations of p 0 . In that case, dw/dx ^ 0 due 
to a gradient dpjdx ^ 0 but du/dx - 0. Physically it can be 
reasoned that u is a better measure of MC than w when it 
comes to bound water diffusion because water molecules 


which is a modified Fick’s second law. The boundary condi¬ 
tions are specified values of u or g or convective conditions 
(Eq. 3). The initial condition is u- w 0 (x, y, t 0 ). 


Materials and methods 

A wood sample made of clear sapwood of Norway spruce 
(Picea abies) with the green dimensions of 31, 42, and 
205 mm in the x (radial), y (tangential), and z (longitudinal) 
directions, respectively, was dried in this study. Five sur¬ 
faces of the sample were coated using polyurethane glue 
(Cascol 1809, Casco) and aluminum foil (Fig. 1). The coated 
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surfaces were also thermally insulated using Styrofoam. 
Three temperature sensors were placed within the sample 
at 1,11, and 21 mm depths from the surface, but in different 
positions in the longitudinal direction. One additional sen¬ 
sor was placed on the surface. During drying, the humidity, 
temperature, and speed of the circulating air was approxi¬ 
mately constant at 43% relative humidity, 50°C, and 4 m/s, 
respectively. A Siemens Somatom AR.T. X-ray CT scanner 
was used to capture a density image in the tangential-radial 
plane of the interior at a constant longitudinal position 
every 10 min during drying. 


where u k wi is evaluated by a parabolic extrapolation to the 
surface using the three u { values that are closest to the 
evaporation surface. 

Evaluation of D from the CT data, 2D method 

Mean u(x , y, t) values for cells, 4.1 x 4.1 x 5 mm 3 , in 3D space 
are evaluated from the CT data. The total mass flux, g k mf , 
transferred from the sample at time t = f is calculated by 
using the total mass decrease, the time step At k = P +i - P~ l 
and the area of the convective surface, A surf . The mass flux is 


Evaluation of D from the CT data, ID method 


u(x , t) values for small volumes (0.14 x 0.14 x 5 mm 3 ) in 3D 
space (voxel values) were measured with CT (see Danvind 
and Moren 6 for a description of the method). In order to 
reduce spread and increase accuracy, mean MCs u(x h 4) for 
seven discrete volumes along the x axis with center posi¬ 
tions at x t = 2.0, 5.9, 9.8,13.7,17.6, 21.6, and 25.5 mm for i = 
1 to 7 were evaluated. A value of u at x 8 = 29.2 mm was set 
equal to the value of u at x = x 7 . x 0 = 0 was the surface 
position, and x 9 = 31.0 mm denoted the inner, insulated 
boundary, u values were measured at discrete time points, 
4, between t Q = Oh and t 200 = 100h with a time step of 0.5h, 
except between t = 83 h and t- 95 h and between t = 42.5 h 
and t- 44 h due to malfunction of the equipment. Denoting 
d/dx with a prime, using mass conservation and a central 
difference scheme we get the mass flux gradient, 


g' k = -p 0i U k = -Po, 




where the subindex i denotes the position x z - = 2.0, 5.9 ... 
29.2 mm for i = 1 to 8 and superindex k denotes the time step 
4 = 0, 0.5,1.0 ... h for k = 1 to 199. Integration of Eq. 8 gives 
the surface mass flux 


gsurf = g(0,4) = - j g’dx = -^g' k AXi (9) 

0 i =1 

where Ax t = x i+1 - x t is the length in the x direction of the 
volume associated with each value u t . The mass flux at posi¬ 
tion i is 


x 1 

gf = g(xi,t k ) = gfurf + j g' k dx = g k ml +^gi k Ax, (10) 

0 /=1 

where Ax, is the length in the x direction of each of the 
volumes that has x coordinates lower than x,. Now Eq. 1 
with a central difference scheme gives 



The mass transfer coefficient /3 is evaluated from Eq. 3 as 
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where subindices i and j denote the positions in the x and y 
directions, respectively. M-l and N -9 are the number of 
cells in the x and y directions, respectively (Fig. 1), and m is 
mass. 

The mass transfer coefficient, fi k , at time t = f is calcu¬ 
lated as 



where u k mij is the MC at time t = £ for the surface of a cell 
adjacent to the evaporation surface with y = y 7 . u k wij is 
linearly extrapolated from the MCs of the two cells next to 
the surface. A smf j is the evaporation surface area of the cell 
at y = y r The modified Fick’s first law of diffusion (Eq. 4) is 
used for the mass flux of the interior cell’s boundary sur¬ 
faces. In the first iteration of an iterative scheme, the cou¬ 
pling term, D xy in Eq. 5, is set equal to zero. For a cell on the 
evaporation surface at position (/, j) = (surf,/), surface mass 
transfer is assumed to govern the mass flux, which is calcu¬ 
lated using p k from Eq. 14. The modified Fick’s second law 
(Eq. 7) is used for mass balances in cells, and the spatial 
derivatives of D are assumed to be zero in a cell. Based on 
these assumptions, Eq. 7 can be stated in numerical form as: 
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In Eq. 15 it is assumed that all elements have the same size, 
Ax, and Ay. From the CT experiments, all parameters except 
D xij and D vij can be achieved for each cell of 30 x 30 voxels 
in the CT image. Hence, there is one equation (Eq. 15) and 
two unknowns for each cell; thus, the equation system for 
all elements is underdetermined. By setting D xij and D y ij 
equal in adjacent cells, in this case 2x6 cells (8.1 x 24.3 mm), 
the number of unknowns for the equation system is re¬ 
duced, and thereby an overestimated system is achieved 
which is solved in a least-squares sense using standard 
Matlab 10 routines. By moving the position of the 2x6 cells 
and repeating the calculation procedure 12 times, once for 
each movement of the cells, 12 sets of D xij and D XII are 
achieved in each cell. The median values of the 12 solution 
sets are then taken as estimated solutions of D x ,, and D v ,, in 
each cell. 

In a second iteration using u Uj , D xij , and D yiJ , the mass 
fluxes, g xii and g yij , are calculated using a numerical form of 
Eq. 4. g xij and g yij are used to recalculate the diffusion 
coefficients, including the coupling term, D xy . This proce¬ 
dure gives two equations and three unknowns, D xx , D yy , and 
D xy , per cell; i.e., when setting up a matrix system of equa¬ 
tions for all cells, an underdetermined system is obtained. 
This time the number of unknowns is reduced by setting 
D xx , D yy , and D xy constant in 2 x 2 cells and by solving the 
overdetermined system in a least-squares sense. As before, 
by moving the position of the 2x2 cells and by solving the 
system four times, four sets of diffusion coefficients per cell 
are achieved in each cell, and the median values are chosen 
as solutions. Using the coupled diffusion coefficients, local 
mass fluxes can be derived, and the procedure of deriving 
diffusion coefficients can be iterated until values stabilized. 
Here, only the first and second iterations have been done. 
In this study, the diffusion in the radial direction (x) was to 
be studied. Therefore, diffusion coefficients in the horizon¬ 
tal direction (y) of the CT image are not presented. 

FEM calculation 

ABAQUS was used with a 3D isothermal and isotropic 
diffusion model as a solver of the diffusion equation (Eq. 6). 
As boundary conditions, we had zero mass flux on five 
surfaces and convection on the evaporation surface (Fig. 1). 
The initial condition u(x, t 0 ) at t Q = 30.1 h was taken from the 
measured data. Based on experimental findings shown be¬ 
low, two alternatives for D were used. As the first alterna¬ 
tive D -D{u) was adapted, and as a second alternative D = 
D(u, x) was used. The measured values of (3(u ) were used 
during the FEM calculations. 


Results 

All results shown in this article relate to the diffusion part of 
the total drying process and hence use a time scale that 
starts at time 30.1 h. The temperature differences between 
three internal positions in the test sample and in the air 



x (mm) 

Fig. 2. Measured basic density p 0 {x) 



Fig. 3. Measured moisture content u(t) for different x values in con¬ 
secutive order from the lowest curve with the lowest x value 


were <0.8°C during the test. p 0 (x) is shown in Fig. 2, u(t ) for 
different x values is shown in Fig. 3, and u{x) for different 
times in Fig. 4. D(u ) for the ID method is shown in Fig. 5 for 
x < 13.7 mm and (3(u ) is shown in Fig. 6. For the 2D method, 
D(u ) is shown in Fig. 7 and (3(u ) in Fig. 6. The objective in 
the FEM calculations was to find D values that gave u(x) 
good correlation to the experimental result at the final time 
t = 99 h. The agreement between experiment and calculation 
was then checked at an intermediate time t - 65 h. The FEM 
simulation using the first alternative D = D(u) according to 
Fig. 8a gives an agreement with the CT measurements 
according to Fig. 4a. The second alternative D = D(u,x) 
according to Fig. 8b has an agreement according to Fig. 4b. 
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a x (mm) 



b x (mm) 


Fig. 4a,b. Moisture content u(x) at different times t measured with 
computed tomography and calculated with finite element method 
(FEM). a FEM approach D = D(u). b FEM approach D = D(u , x) 




Fig. 5a,b. D(u ) evaluated with the one-dimensional (ID) method for 
different x values, ar = 4.0 and 5.9mm. b x = 9.8 and 13.7mm 


The initial conditions u 0 at t Q = 30.1 h for the FEM simula¬ 
tions are the same for both alternatives. 


Discussion 

For the ID method, D shows a rather clear dependence on 
u and x as is seen from Fig. 5. The curves for x = 5.9, 9.8, and 
13.7 mm agree well (Fig. 5) but the curve for x = 4.0 mm is 
translated to the left and downward compared with the 
other curves. The curve for x = 17.6 mm (not shown) agrees 
quite well with the curves for x = 5.9, 9.8, and 13.7 mm, but 
shows more spread. The curves for x = 21.6 and 25.5 mm 
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(not shown) show even more spread, but they seem to agree 
with the curves for x = 5.9, 9.8, 13.7, and 17.6mm. Thus, 
there is a unique curve for x = 4.0 mm, and all of the rest of 
the curves at the other x values seem to agree reasonably 
well. The reason for the large spread of D for large x values 
is probably the decreasing MC gradient and the decreasing 
mass flux with depth (x), which make the numerical errors 
large (Eq. 11). The overall conclusion is that the diffusion 
coefficient is a function of MC and depth, D = D(u,x ). The 
dependence on u is especially large in the interval where 
16% < u < 30%, while for u < 16%, D seems rather constant. 


-4 

x 10 



Fig. 6. p(u) evaluated with the ID and 2D methods with alternative 
surface moisture contents. w ext , linear extrapolation from two adjacent 
u values or parabolic extrapolation from three adjacent u values; u mid , 
mean value of the surface cell; i.e., u mid = u(x i= i, r) 


The dependence on the distance from surface x is only 
significant when x < 5mm, i.e., near the surface. p shows 
more spread at lower mean MC due to smaller w surf - than 
at higher mean MC (Fig. 6). 

For the 2D method, only D x values are presented. This is 
due to large spread in D y , probably caused by the one- 
directional drying resulting in almost constant values of u in 
the y direction, which strongly influenced the derivation of 
D y . D x values (Fig. 7) are similar to the results from the ID 
method, but the results have less spread. A reason for the 
lower spread is probably the calculation method in which 
the median values of several solution sets are taken as D. 
P values for the 2D method are lower than for the ID 
method (Fig. 6). This is due to the different choices of u smf . 

When trying to reproduce the original, measured u(x ) 
values at t - 99.1 h with the FEM calculation, the second 
alternative with D = D(u , x) is best (Fig. 4b). This shows the 
validity of the experimentally derived D values and the 
dependence of D on distance to the evaporation surface. 
The first alternative with D = D(u) gave good correlation at 
t = 99.1 h but poor correlation at t- 65.1 h (Fig. 4a). 

In Fig. 8b, evaluations of D from Hukka 1 and Rosenkilde 
and Arfvidson 4 are compared with our values used for FEM 
calculations. Their values and ours agree quite well, at least 
for u < 15%. The discrepancy in D for u > 20% between our 
values and theirs could be due to differences in wood mate¬ 
rial (Hukka 1 used Norway spruce heartwood and 
Rosenkilde and Arfwidson 4 used Scots pine sapwood) and 
differences in evaluation methods. Hukka 1 assumed that 
D(u) is an exponential function and Rosenkilde and 
Arfvidson 4 used another type of curve-fitting method. The 
rate of the mean MC for a wood sample is essentially con¬ 
trolled by P and not D when u is high, and vice versa. 
However, an appropriate D(u ) description is important for 
realistic local u(x) values. 

The dependence on depth can be a dependence not on 
depth itself but via some other parameter (not measured 


Fig. 7. Diffusion coefficient 
D(u,x ) evaluated with the 2D 
method, x = 2.0, 5.9, 9.8,13.7, 
17.6, 21.6, 25.5mm shown as 
alternating series of dots and 
crosses 
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II (%) 



b u (%) 

Fig. 8. a D(u) used in the first approach of FEM calculations, b D(u, x ) 
used in the second approach of FEM calculations. Linear interpolation 
of D for 0 < x < 8 mm using the curves for x = 0 and x = 8 mm. 
Comparison values of Hukkas 1 values are for Norway spruce heart- 
wood, those of Rosenkilde and Arfvidssons 4 are for Scots pine 
sapwood 


here), which itself is a function of depth. Such a parameter 
is stress, because there are probably large stresses near the 
surface. Also, u is a function of depth because the material 
near the surface has a faster changing moisture situation 
than the material far away from the surface. 

Another possible reason for the behavior of D near the 
surface could be the position of the evaporation front dur¬ 
ing drying. Wiberg 7 and Rosenkilde and Arfvidsson 4 show 
how the evaporation front recedes into the material in the 
capillary regime of drying, creating a dry “shell” near the 
surface. This dry shell probably exists also in the beginning 
of the diffusion regime. The dry shell may have caused u i=l j 
to have lower values than they would have had if the dry 


shell had not existed. Hence, one can propose to include a 
dry shell formulation in future evaluations of D, as sug¬ 
gested by Salin. 12 

The calculations of D are sensitive to measurement er¬ 
rors of p 0 , u , t , and x values, because derivatives of u in space 
and time are estimated by numerical schemes (Eqs. 8-11 
and Eqs. 13-15). An error estimate of u based on the spread 
of the graph of ii(t) calculated with a central difference 
approximation (Eq. 8) was made. The assumption was that 
u(t ) is a normally distributed stochastic variable and that t 
values are exact, which results in an estimated standard 
deviation of u of the order of 0.04% when u is approxi¬ 
mately 20%. This standard deviation of u is considered low 
in comparison with earlier values (see Danvind 13 ). It is 
believed that this spread in measured u values is the main 
cause of spread in D. The spread in D increases when u 
decreases. 
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